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Abstract 

The control of complex systems and network-coupled dynamical systems is a topic of vital 
theoretical importance in mathematics and physics with a wide range of applications in engi¬ 
neering and various other sciences. Motivated by recent research into smart grid technologies 
we study here control of synchronization and consider the important case of networks of cou¬ 
pled phase oscillators with nonlinear interactions-a paradigmatic example that has guided our 
understanding of self-organization for decades. We develop a method for control based on iden¬ 
tifying and stabilizing problematic oscillators, resulting in a stable spectrum of eigenvalues, and 
in turn a linearly stable synchronized state. Interestingly, the amount of control, i.e., number 
of oscillators, required to stabilize the network is primarily dictated by the coupling strength, 
dynamical heterogeneity, and mean degree of the network, and depends little on the structural 
heterogeneity of the network itself. 


Introduction 

Complex networks and complex systems describe the physical, biological, and social structures that 
connect our world and host the dynamical processes vital to our lives HIM. The failure of such 
large-scale systems to operate in the desired way can thus lead to catastrophic events such as power 
outages M , extinctions M, and economic collapses M- Thus, the development and design 
of efficient and effective control mechanisms for such systems is not only a question of theoretical 
interest to mathematicians, but has a wide range of important applications in physics, chemistry, 
biology, engineering, and the social sciences mm- 

The roots of modern linear and nonlinear control reach back several decades, but recently re¬ 
search in this direction has seen a revival in physics and engineering communities. For instance, 
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the concept of structural controllability , which is based on the paradigm of linear homogenous dy¬ 
namical systems, was first introduced by Lin in [T2] and more recently investigated in mm- 
These advances have enabled further progress related to structural controllability such as cen¬ 
trality m, energy [IB], effect of correlations m emergence of bimodality [IB], transtion and 
nonlocality [19] , the specific role of individual nodes m , target control [21], and control of edges 
in switchboard dynamics [ 223 - Significant advances have also been made in the control of nonlinear 
systems, for instance the control of chaotic systems using unstable periodic orbits [23] . control via 
pinning [23] [253, [26] , control and rescue of networks using compensatory perturbations mm, 
and control via structural adaptation [29] . Implicit in all such network control problems are the 
questions of (i) what form(s) of control should one choose? and (ii) how much effort is needed to 
attain a desired state [30]? 

Motivated by ongoing studies on the stability and function of power grids ED E2], we study 
here the control of heterogeneous coupled oscillator networks [331 EIJ . Recent research into smart 
grid technologies has shown that certain power grid networks called microgrids evolve and can 
be treated as networks of Kuramoto phase oscillators [35] . A microgrid consists of a a relatively 
small number of localized sources and loads that, while typically operating in connection to a 
larger central power grid, can disconnect itself and operate autonomously as may be necessitated 
by physical or economical constraints. In particular, by means of a method known as frequency- 
drooping , the dynamics of microgrids become equivalent to Kuramoto oscillator networks - a class 
of system for which a large body of literature detailing various dynamical phenomena exists [36], 
Here we develop a control mechanism for such coupled oscillator networks, thus providing a solution 
with potentially direct application to the control of certain power grids. 

Our goal is to induce a synchronized state, a.k.a. consensus, in a given coupled oscillator net¬ 
work and guarantee asymptotic stability by applying as few control gains to the network as possible. 
Our method is based on calculating the Jacobian of the desired synchronized state and studying its 
spectrum, by which we identify the oscillators in the network that contribute to unstable eigenval¬ 
ues and thus destabilize the synchronized state. Importantly, our method not only identifies which 
oscillators require control, but also the required strength of each control gain. Interestingly, we 
find that the control required to stabilize a network is dictated by the coupling strength, dynamical 
heterogeneity, and mean degree of the network, and depends little on the structural heterogeneity 
of the network. In other words, the number of nodes requiring control depends surprisingly little 
on the network topology and degree distribution and is more sensitive to the average connectivity 
of the network and the dynamical parameters. While Kuramoto oscillator networks serve as our 
primary system of interest due both to its specific correlation with mircogrids as well as its rich 
body of literature, we note that our method can be applied to a much wider set of oscillator net¬ 
works, provided that their linearized dynamics take a certain form. Moreover, since Kuramoto and 
other oscillator network models have served as a paradigmatic example for modeling and studying 
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synchronization in various contexts, we hypothesize that our results may shed light more generally 
on the control of synchronization processes and could potentially give insight into other important 
applications such as the termination of cardiacarrhythmias m and treatments for pathological 
brain dynamics [5H J. 

Results 

The Kuramoto model 

We consider the famous Kuramoto model for the entrainment of many coupled dissipative oscil¬ 
lators [32]. The Kuramoto model consists of N phase oscillators 0* for i = 1 that, when 

placed on a network dictating their pair-wise interactions, evolve according to 

N 

9i = Ui + Aij sin (Oj - 0j). (1) 

j = 1 

Each oscillator i has a unique nature frequency that describes its preferred angular velocity in the 
absence of interactions, which is typically drawn randomly from a distribution g(u). Furthermore, 
the global coupling strength K describes the influence that oscillators have on one another via 
the network connectivity, which is encoded in the adjacency matrix [d,j]. Here we focus on the 
simple case of an undirected, unweighted network (Ajj = 1 if oscillators i and j are connected 
by a link and Aij = 0 otherwise), but we note that all results presented here easily generalize to 
directed and weighted networks. We also assume that the network is connected, i.e., irreducible. 
Over the last few decades the Kuramoto model has proven to be very useful for modeling real- 
world systems [36] 0Q], uncovering the mechanisms behind emergent collective behavior EH 02 , 
exploring additional effects such as time-delays [33] and community structure [33J, and finding 
optimal networks structure [45] , 

Depending on the coupling strength K, as well as the frequency vector u> and the network 
topology, the steady-state dynamics of Equation (JT]) can attain many different states that included 
complete incoherence, partial synchronization, and full synchronization. The latter is characterized 
by lim^oo \9j(t) — 9{(t)\ = 0 and is also referred to as full phase-locking, frequency-synchronization, 
or consensus. The fully synchronized state (henceforth called the synchronized state) typically 
displays a large degree of phase synchronization r ~ 1, where re= iW 1 YljLi e ldj is the standard 
Kuramoto order parameter. In Figure Oja) we illustrate a synchronized state in a group of five 
oscillators, each moving with an angular velocity of w. The order parameter is illustrated as vector 
of length r with an offset angle ^ from the positive real axis. 
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Figure 1: Control of coupled oscillator networks. For a five-oscillator 
network, illustration of (a) the fully synchronized state with the Kuramoto 
order parameter, (b) the Jacobian of a stable synchronized state, and (c) control 
applied to oscillators i = 2 and 3. 


Stability and control 

In this Article we address the problem of control by first assuming that due to the system parameters 
(i.e., the coupling strength, natural frequency sequence, or network topology), the steady-state 
dynamics of Equation (JT]) are at least partially incoherent, i.e., one or more oscillators remains 
desynchronized. In the example of the power grid, a single desynchronized oscillator represents 
a single power failure, but can have further damaging effects, in particular triggering a cascade 
of additional failures and ultimately a power outage [5j. Thus, our goal is to find a synchronized 
state and stabilize it. If all oscillators are initially synchronized, then our goal is trivially realized, 
however our method can be used to make this state more robustly stable. In a synchronized state 
like the one we seek here we expect the oscillators to be clustered in a single reasonably tight cluster 
such that |Oj — 9{\ <C 1, and thus Equation P) can be linearized to 

N 

t (2) 

j =i 

where L is the network Laplacian matrix with entries = 5ij ^ An — Aij and Sij is the Kronecker 
delta. A straight-forward analysis yields a “target” synchronized state (within the rotating reference 
frame 9 H > 9 + ( uj)t ) given by the vector 6* = K~ 1 L^UlJ, where L t is the pseudo-inverse of the 
Laplacian [35]. (We summarize a derivation of this result in Materials and Methods.) We note 
that, since the system is assumed to be partially incoherent, the fixed point 6 = 6* either does not 
exist or is unstable. However, we take 6 = 6* to represent the closest synchronized fixed point for 
the given parameter values, and therefore we use it a target. We also note that, while Equation ([2]) 
was directly obtained from linearizing Equation P), other systems of more general forms yield 
equivalent linearizations and therefore can also be controlled using the method we provide here. 
We present an example of such a general system with arbitrary coupling function (e.g., see m) in 
the Materials and Methods. 

The stability of 6 = 6* is dictated by the Jacobian matrix whose entries are defined [DF]ij = 


4 






















ddi/ddj, and is stable if all the eigenvalues of DF\g* are non-positive. In our case we have that 


DFij — 


~ K Hi# A a cos(6»* - Of) 
KAij cos {0* - Of) 


if i=j 

otherwise. 


( 3 ) 


We note that each row (and column) of DF sums to zero, i.e., satisfies DFa = — DFij. This is 
a particularly convenient property for using the Gershgorin circle theorem [38], which implies that 
eigenvalues of DF lie within the union of closed discs Di for i = 1,..., N, which are each centered 
at DFa and have radius Ri, where Ri = \DF t] \. (The full theorem is given in Materials and 

Methods.) In particular, if all the off-diagonal entries of DF are non-negative, then it follows that 
each Gershgorin disc is contained in the left-half plane, implying that all eigenvalues are non-positive 
and the solution is stable. An illustration of this case is presented in Figure [2(b). If, however, one 
or more non-diagonal entries of DF are negative, then each Gershgorin disc corresponding to a row 
with a negative off-diagonal entry enters the right-half plane, admitting the possibility for one or 
more positive eigenvalues and thus destabilization. Thus, the oscillators that require control can be 
easily identified as those whose corresponding rows have one or more negative off-diagonal entries. 

We aim to stabilize the synchronized solution by adding one or more control gains to the system, 
as illustrated in Figure [Die). In following recent literature, we will refer to oscillators to which we 
apply control as driver nodes, and to oscillators to which we do not apply control as free nodes. 
We choose the control gains to take the form fi(t) = F)sin(0j — Of), where F) is the strength of 
the i th control gain and fa is a target phase that can in principle depend on either local or global 
information, and vary in time. Here we focus on the choice of target phase fa = 6*, and discuss 
other possibilities below. Since the control gain depends on the current state of the system, this 
can be though of as a form of feedback control. The new dynamics are then given by 


N 

9i = uh + A ij sin (Oj ~ Of) + F 1 * sin(0* - Of), (4) 

3 = 1 

where we take F) = 0 for free nodes. While the off-diagonal entries of DF remain unaltered, the 
new diagonal entries are given by DFa = —K Aij cos (Oj —0*) — F). Thus, we set coupling gain 
strength of each driver node i such that it satisfies F) > K Aij [| cos (0* — Of) | — cos (0* — Of)]. 
This ensures that all Gershgorin discs are contained in the left-half plane, implying that (up to the 
linear approximation of 0*) all eigenvalues are non-positive and the synchronized state is stable. 
(In the case of a directed network this implies that all eigenvalues have non-positive real part and 
the synchronized state is stable.) 

We now briefly comment on the choice of target phases fa in the control gains. In the method 
outlined above, we have set the target phase equal to the steady-state phase, fa = Of. This is a 
convenient choice for the derivation described above. Additionally, we find that in practice other 
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Figure 2: Control of coupled oscillator networks: random networks. 

Example of control applied to a coupled oscillator network with ER and SF 
topologies (top and bottom rows, respectively) with (k) = 6 and N = 1000. 
The coupling strength is K = 0.4. Time series for phases 9i(t) of 10% of 
oscillators (a), (d) without and (b), (e) with control. Driver nodes constituted 
37.1% and 44.5% of the ER and SF networks, respectively, (c), (f) Time series 
for the degree of phase synchronization r(t) for the networks without (red) and 
with (blue) control. 


choices also yield positive results. In particular, one choice that tends to yield slightly better results 
is to force each driver node towards the center of the synchronized cluster, i.e., 4>i = cj), where we 
assume the cluster is centered at the angle <fi. Target phases can also be chosen according to global 
or distributed control strategies. In particular, given the standard Kuramoto order parameter 
re^ ■ e i9 i /N or the set of local order parameters r l e i ^ i = Aije ldj , the choices <pi = ip 

and (j)i = ipi correspond to global and distributed control strategies that typically yield favorable 
results. We also note before presenting examples that since the method outlined above depends on 
the approximation of the steady-state solution 6* ~ in practice we add a buffer margin 

when identifying unstable oscillators, looking for non-diagonal entries of DF that are not necessarily 
negative, but rather less than some e > 0. We find that the choice e = 0.2 is sufficient, and is what 
we use in the examples below. 

Control of random networks 

We now demonstrate our approach by considering two types of random networks: Erdos-Renyi [39] 
(ER) networks and scale-free (SF) networks. Each ER network is constructed using a fixed link 
probability p and each SF network is built using the configuration model m with a degree sequence 
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drawn from the distribution P(k ) oc k ~ 7 with 7 = 3 and enforced minimum degree ko- To tune 
the mean degree (k) of each network we set either p = (k)/(N — 1) or ko = (k)/( 7 — 1). Figure [2] 
illustrates our results with an example of each type of network, where we have used networks of 
size N = 1000 with mean degree (k) = 6, set the coupling strength K = 0.4, and used natural 
frequencies drawn from a uniform distribution with zero mean and unit variance. The top and 
bottom rows display the results for the ER and SF networks, respectively, displaying the time 
series 9{t) of a randomly selected 10% of the oscillators without control [first column, panels (a) 
and (d)] and with control [middle column, panels (b) and (e)], after discarding a long transient. 
The difference between no control to control is quite drastic, with a large fraction of desynchronized 
oscillators without control, and full synchronization with control. Driver nodes constituted 37.1% 
and 44.5% of the ER and SF networks, respectively, to attain the synchronized state. In the last 
column [panels (c) and (f)] we present the degree of phase synchronization, plotting the order 
parameter r{t) for both solutions with and without control. Unlike the solution without control, 
which fluctuates significantly at a relative low value, the solution with control reaches a steady 
value near r[t) ~ 1. 

Next we investigate the properties of driver and free nodes by revisiting our method. This is 
an essential question due to the heterogeneity of the oscillators both in terms of network structure 
(i.e., degree distribution) and also local dynamics (i.e., natural frequencies). An oscillator i is a 
driver node if one of its off-diagonal entries DF t] oc cos(6j — 9*) is negative, and therefore oscillators 
with large (small) steady-state values 6* gc |[LW|j| tend to be driver (free) nodes. Furthermore, 
we find that these values scale approximately linearly with the ratio of the natural frequencies to 
degrees, i.e., [Ltuj\i cu i/ki. We illustrate this is Figured where we plot the relationship [L^u\i 
vs uJi/ki for the example ER and SF networks presented above [panels (a) and (b), respectively], 
denoting driver nodes with red crosses and free nodes with blue circles. These results show the 
important role that dynamics, in addition to network structure, plays in dictating controlling the 
system. In particular, driver nodes of the system tend to balance a large ratio (in absolute value) 
of natural frequencies to degrees. 

Finally, we quantify the overall effort required for consensus by studying how the fraction of 
driver nodes, denoted nr> = Np/N, where Njj is the total number of driver nodes, depends on 
both the system’s dynamical and structural parameters. Presenting our results in Figure 01 we 
first explore how the fraction of driver nodes depends on the coupling strength by plotting in panel 
(a) no vs K for both ER and SF networks with mean degrees (k) = 4, 8, and 12 (blue circles, 
red triangles, and green squares, respectively). Results for ER and SF networks are plotted with 
unfilled and filled symbols, respectively, and each curve represents an average over 100 network 
realizations, each averaged over 100 random natural frequency realizations. While it is expected 
that no decreases monotonically with K, the curves’ dependence on network topology and mean 
degree is nontrivial. In particular, the the shape of no vs K depends more sensitively on the mean 
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Figure 3: Driver and free nodes. \L'u\i vs u)i/ki for driver (red crosses) and 
free (blue circles) nodes in the (a) ER and (b) SF networks used in Figure [2] 


degree than the topology, suggesting that network heterogeneity has little effect on overall control 
in comparison to average connectivity. In light of the significant dependence of overall control on 
the coupling strength, we investigate the coupling strength required to synchronize a network if 
limited mount of control is available. To this end we calculate for each family of networks the 
required coupling strengths K 5 %, K w %, and K 2 q% for which, on average, a fraction no = 0.05, 0.1, 
and 0.2 will achieve synchronization as a function of the average degree (k). We plot the results in 
Figured] (b). We point out again that ER and SF networks behave very similarly on average, and 
that with a larger mean degree, a smaller coupling strength is required to achieve synchronization. 


Discussion 

Theoretical and practical aspects of the control of dynamical processes remains an important and 
ongoing area of interdisciplinary research at the intersection between mathematics, physics, biology, 
chemistry, engineering, and the social sciences. Control of complex networks and complex systems 
is particularly important since together they comprise most of the world we live in eh, however the 
nonlinear nature of realistic dynamical processes and the complex network topologies of real net¬ 
works represent challenges for the scientific community. Building on concepts from classical linear 
control theory ||2], recent work has made significant advances in understanding structural control¬ 
lability Banai, and significant progress has been made in the development of control mechanisms 
for networks of nonlinear systems [23j i2H ESj- Nonetheless, due to the problem-sensitive nature of 
most real-world problems and applications requiring control techniques, further progress in design¬ 
ing and implementing efficient and effective control mechanisms for a wide range of problems with 
practical constraints remains an important avenue of research. 

In this Article we have focused on the control of synchronization, i.e., consensus, in coupled 
oscillator networks. Our primary inspiration has been advances in the research of power grid 






Figure 4: Fraction of driver nodes, (a) The fraction of driver nodes no in 
ER (unfilled symbols) and SF (filled symbols) as a function of coupling strength 
K for mean degrees (k) = 4, 8, and 12 (blue circles, red triangle, and green 
squares, respectively), (b) The required coupling strengths K 5 %, K iq%, and 
K 2 q% required to achieve consensus given no = 0.05, 0.1, and 0.2 (blue circles, 
red triangle, and green squares, respectively). Each data point is the average 
over 100 network realizations of size N = 1000, each averaged over 100 natural 
frequency realizations. 


networks [52, [53]. In particular, recent studies have shown that certain power grids known as 
microgrids can be treated as Kuramoto oscillator networks [351 [351 . Here we have presented a control 
method that can easily be applied to Kuramoto networks and other phase oscillator networks, thus 
providing a control framework with potentially direct application to these new technologies. Our 
method is based on identifying and stabilizing a synchronized state for a given network via spectral 
properties of the Jacobian matrix and we have demonstrated its effectiveness on both Erdos-Renyi 
and scale-free networks. We have observed that driver nodes, i.e., oscillators that require control, 
tend to balance (in absolute value) large natural frequencies with small degrees. Furthermore, the 
overall amount of control required to achieve synchronization decreases with both coupling strength 
and mean degree, while the total effort required to attain a synchronized state depends sensitively 
on the average connectivity of the network and the dynamical parameters, but surprisingly little 
on the network topology and degree distribution. These results enhance our understanding of and 
ability to understand, optimize, and ultimately control synchronization in power-grid networks (see 
in particular um), and more generally complement important work on the control of network- 
coupled nonlinear dynamical systems [261 f28l .29]. 

While our central inspiration and target application is in the area of power grid technology, 
synchronization phenomena plays a vital roll in a variety of complex processes that occur in both 
natural and man-made systems, including healthy cardiac behavior [53], functionality of cell cir¬ 
cuits [55], stability of pedestrian bridges [55], and communications security [57]. Given this broad 
range of applications, we hypothesize that our findings here may potentially shed some light on 
the control of synchronization in other contexts, for instance cardiac physiology and neuroscience. 
For instance, a large amount of research has recently been devoted to the development of car- 
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diac arrhythmia treatments that require minimal shock to knock out fatal asynchronous behavior 
such as cardiac fibrillation [58| and the promotion of normal brain oscillations [59] while repressing 
disorders such as Parkinson’s disease which are associated with abnormal oscillations [ 60] . 


Materials and Methods 

Steady state solution 

To derive the steady state solution 6* = K~ l L^u, we begin with Equation Q, which represents the 
linearized dynamics of Equation (JT]). Recall that this linearization requires that we are searching 
for a synchronized state where all oscillators are tightly packed in a single cluster so we expect that 
1 6j —0i\ ^ 1- We also note that the mean frequency of all oscillators is given by the mean natural 
frequency (w). For simplicity we enter the rotating frame 6^-0 + (c u)t, effectively setting the mean 
frequency to zero. It is then convenient to write Equation ([2]) in vector form, i.e., 

KLO , (5) 

where L is the network Laplacian whose entries are defined Ljj = 5ij ^ An—Aij. While L has a zero 
eigenvalue, denoted Ai = 0, rendering it non-invertible, it does have a pseudo-inverse defined using 
its other eigenvalues (which are non-zero provided that the network is connect) and corresponding 
eigenvectors, L 1 = |46j . Each eigenvector is normalized such that { v i7 }jL 2 f° rms 

an orthonormal basis for the space of vectors in with zero mean. Thus, both L and L 1 share 
a nullspace which is spanned by the eigenvector v 1 oc 1, and therefore map vectors onto the space 
of zero-mean vectors in R w . With the pseudoinverse in hand, we can finally obtain the desired 
steady-state solution by setting 0 = 0 and solving for 0, which yields the solution 6* = K 
as desired. 

General oscillator networks 

Here we present an example of a more general oscillator network than that in Equation (]T]) that 
can be controlled using the same method detailed above. In particular, we generalize to account 
for an arbitrary coupling function H(0), yielding 

N 

0 i = u i + Kj2 AijH(0j - 0i). (6) 

3 = 1 

We assume that H{0) is 27r-periodic and at least once continuously differentiable. Importantly, H 
need not satisfy H{ 0) = 0, and thus coupling between neighboring oscillators can be frustrated [47], 
denoting that even when two oscillators are exactly equal, their interaction term does not vanish. 
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Provided that the coupling frustration is not too large, e.g., H(0)/V2H'(0) <l,a tightly clustered 
synchronized state is attainable, and linearizing Equation (J6|) yields 

N 

9i^uJi + KH(0)h - KH\ 0) ^ Lijdj. (7) 

3 = 1 

Importantly, by dehning the quantities Ui = uii + KH{0)kj. and K = KH'( 0), it is easy to see that 
the linearized dynamics of Equation ([7]) are of the same form as Equation (|2|), and therefore the 
control method we present above can be readily applied. 

Gershgorin circle theorem 

Definition. (Gershgorin Discs) Let M be an N x N complex matrix. For i = 1,..., N let Ri = 
-- \Mij\ be the sum of absolute values of non-diagonal elements of row i, and define D(Ma,Ri ) 
closed disc of radius R{ centered at Mu. Di = D(Mu,Ri ) is the i th Gershgorin Disc. 

Theorem. (Gershgorin) All eigenvalues of the matrix M lie within the union U f = \Di of Gershgorin 
discs. 
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